snakemake 使用多环境分析数据

原创 JunJunLab 老俊俊的生信笔记 2022-08-15

收录于合集

#100天生信/数据科学自我挑战 100 个

#RNA-seq 39 个

看万山红遍,层林尽染

1引言

前面已经有两期内容讲了 snakemake 的基本用法,基本上我们 写个简单的 workflow 是完全没有问题的。

我们知道当我们分析数据走一个完整的 pipeline 时,可能每个 step 用到的 软件不一样 ,每个软件的 依赖的软件也不一样 ,此外各种软件的 依赖环境也可能不一样 。

举个例子: 有些软件只能在 python2 的环境运行,到了 python3 就会报错。

一个很常见的现象就是你在一个环境里安装了太多软件,可能会导致之前的软件无法运行或者新装的软件也可能会报错。

最好的办法就是建立不同的环境来安装这些软件方便管理。 snakemake 可以使用 conda 参数 来提供你的 环境 yaml 文件 ,从而达到每个模块可以在不同的环境进行分析。

2创建 yaml 文件

假如我们创建一个数据预处理的 python3 环境:

$ conda create -n prefilter python=3 fastqc multiqc fastp

激活环境:

$ conda activate prefilter

导出环境:

$ conda env export > prefilter.yaml

# 查看内容
$ head -n 15 prefilter.yaml
name: prefilter
channels:
  - http://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda
  - http://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/main
  - http://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/main/
  - http://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/
  - http://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda/
  - defaults
dependencies:
  - _libgcc_mutex=0.1=main
  - _openmp_mutex=4.5=1_gnu
  - blas=1.0=mkl
  - bzip2=1.0.8=h7b6447c_0
  - c-ares=1.17.1=h27cfd23_0
  - ca-certificates=2021.10.26=h06a4308_2
...

这个 yaml 文件包含了软件及依赖软件的具体信息,还有软件下载来源。

3使用 yaml 文件

那么如何在 snakemake 的 workflow 中使用呢? 我们需要在模块里引入 conda 参数,写一个 fastqc 和 multiqc 的小模块:

rule filter_qc:
    input:
        expand("fastq-data/{sample}.fq.gz", sample=samples)
    output:
        expand("QC-data/{sample}_fastqc.html",sample = SAMPLES),
        expand("QC-data/{sample}_fastqc.zip",sample = SAMPLES)
    conda:
        "envs/prefilter.yaml"
    shell:
        "fastqc -t 10 -o QC-data/ {input}"

rule multiQC:
    input:
        expand("0.QC-data/{sample}_fastqc.zip",sample = SAMPLES)
    output:
        "0.QC-data/multiqc_report.html"
    conda:
        "envs/prefilter.yaml"
    shell:
        "multiqc -o 0.QC-data/ {input}"

在执行时，Snakemake 将 自动创建该环境，并在其中执行 shell 命令,我们讲 yaml 文件放在了 envs 文件夹下。

测试:

$ snakemake --use-conda -n

4官网示例代码

最后提供一下官网的示例完整 pipeline 代码:

samples = ["A", "B"]


rule all:
    input:
        "calls/all.vcf",
        "plots/quals.svg"


rule bwa:
    input:
        "data/genome.fa",
        "data/samples/{sample}.fastq"
    output:
        temp("mapped/{sample}.bam")
    conda:
        "envs/mapping.yaml" ## 比对环境
    threads: 8
    shell:
        "bwa mem -t {threads} {input} | samtools view -Sb - > {output}"


rule sort:
    input:
        "mapped/{sample}.bam"
    output:
        "mapped/{sample}.sorted.bam"
    conda:
        "envs/mapping.yaml" ## 比对环境
    shell:
        "samtools sort -o {output} {input}"



rule call:
    input:
        fa="data/genome.fa",
        bam=expand("mapped/{sample}.sorted.bam", sample=samples)
    output:
        "calls/all.vcf"
    conda:
        "envs/calling.yaml" ## calling 环境
    shell:
        "samtools mpileup -g -f {input.fa} {input.bam} | "
        "bcftools call -mv - > {output}"

rule stats:
    input:
        "calls/all.vcf"
    output:
        report("plots/quals.svg", caption="report/calling.rst")
    conda:
        "envs/stats.yaml" ## 统计环境
    script:
        "scripts/plot-quals.py"